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Abstract —Transportation processes, which play a prominent 
role in the life and social sciences, are typically described 
hy discrete models on lattices. For studying their dynamics a 
continuous formulation of the problem via partial differential 
equations (PDF) is employed. In this paper we propose a symbolic 
computation approach to derive mean-field PDFs from a lattice- 
based model. We start with the microscopic equations, which 
state the probability to find a particle at a given lattice site. 
Then the PDFs are formally derived by Taylor expansions of the 
probability densities and by passing to an appropriate limit as 
the time steps and the distances between lattice sites tend to zero. 
We present an implementation in a computer algebra system that 
performs this transition for a general class of models. In order 
to rewrite the mean-field PDFs in a conservative formulation, 
we adapt and implement symbolic integration methods that can 
handle unspecified functions in several variables. To illustrate our 
approach, we consider an application in crowd motion analysis 
where the dynamics of bidirectional flows are studied. However, 
the presented approach can be applied to various transportation 
processes of multiple species with variable size in any dimension, 
for example, to confirm several proposed mean-field models for 
cell motility. 

I. Introduction 

Mean-field models play an important role in applied math¬ 
ematics and have become a popular tool to describe transporta¬ 
tion dynamics in the life and social sciences. In the derivation 
of such models the effect of a large number of individuals on a 
single individual is approximated by a single averaging effect, 
the so called mean-field. Applications include cell migration 
at high densities, cf. [1], [2], transport across cell membranes 
as occurring in ion channels, cf. [3], [4], traffic flow [5] as 
well as the motion of large pedestrian crowds, see e.g. [6], 
[7]. Understanding the complex dynamics of large interacting 
groups of particles is of high practical relevance and initiated a 
lot of research in the held of physics, transportation research, 
and applied mathematics. 

Mathematical models on the micro- as well as the macro¬ 
scopic level have been used successfully to describe various 
aspects of these transportation processes. On the microscopic 
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level the dynamics of each individual are modelled taking into 
account its interactions with all others as well as interactions 
with the physical surrounding. This approach results in high¬ 
dimensional and very complex systems of equations. On the 
macroscopic level the crowd is treated as a density which 
evolves according to a partial differential equation (PDE) or 
systems thereof. The transition from the microscopic to the 
corresponding macroscopic description is an active area of 
research with a lot of open analytic questions. 

On the microscopic level a distinction is made between 
two different models: force-based or lattice-based models. In 
the former the dynamics of each individual is determined by 
the forces acting upon it, i.e. exerted from the others and the 
surrounding; the latter states the probability to And a particle 
at a discrete position in space (the lattice point) given the 
transition rates of the particle to move from one discrete lattice 
point to another. 

Lattice-based models, also known as cellular automata, are 
a very prominent tool to describe cell motility, cf. [8], as well 
as pedestrian dynamics (cf. [9], [10]), since exclusion pro¬ 
cesses can be included naturally. In exclusion-based processes 
each lattice site can be occupied by at most one individual, 
giving a simple way to account for the finite particle size. 
In the last years there has been an increasing interest in the 
derivation of the corresponding continuum equations in both 
fields, see for example [1], [2] in case of cell dynamics or [6], 
[7] describing pedestrian dynamics. The general structure of 
the resulting mean-field equations depends on the transition 
rates, but common features include 

1) their conservative nature; i.e. they are based on the 
assumption that the total mass is conserved; 

2) an underlying gradient flow or perturbed gradient 
flow structure with respect to a certain metric; so¬ 
lutions of the first one correspond to minimizers of 
an energy functional with respect to a certain metric. 

These structural features allow to write the mean-field PDEs 
in terms of the diffusivity and energy functionals; quantities 
which are of interest for analysists. Therefore it is desirable 
to derive the mean-field equations in this conservative form, 
although the general formulation is not unique. 

Different strategies have been used to pass to the macro¬ 
scopic limit, i.e. to derive the corresponding continuum equa¬ 
tions as the number of particles tends to infinity, in either 



approach. In force-based models the macroscopic limit can be 
derived using the so called BBGKY hierarchies, see for exam¬ 
ple [ 11 ], initially developed in the field of statistical physics. 
In the case of stochastic underlying dynamics the derivation 
of the mean-field description has been studied rigorously for 
simpler models by considering the hydrodynamic limit, see 
[12]. Simpler models, such as the Patlak-Keller-Segel model 
for chemotaxis or reaction-diffusion equations, were rigorously 
derived for a stochastic many-particle system, see [13] and [14] 
respectively. 

We would like to mention a related work by Penington and co¬ 
workers [15] on a systematic construction method to determine 
the continuum limit of nonlinear PDEs from discrete lattice 
based models. Their approach is based on representing the 
transition rates using appropriate rotation operators as well 
as symmetry conditions to derive general expressions for 
the transportation coefficients in the corresponding nonlin¬ 
ear PDEs. This technique can be used for a large class of 
problems (including multi-species dynamics in various space 
dimension), but assumes that the transition rate of a species 
depends only on the average occupancy of a site by any of the 
different species and not if the site is occupied by a particular 
subpopulation or not. This approach cannot be applied to the 
pedestrian model presented later on, since the transition rates 
depend on the affiliation to either group. 

In the case of a lattice-based model we 

1 ) replace the probability to find a particle at a lattice 
site by a formal Taylor expansion (up to a certain 
order) of the corresponding density, 

2 ) pass to an appropriate limit as the lattice size and 
time tends to zero (dropping higher-order terms). 

In this paper we present an algorithmic approach to 
derive the corresponding continuum equations from a 
lattice-based model using tools from symbolic computation. 
While it is relatively straightforward to perform the formal 
Taylor expansions and the corresponding limit, it is a more 
challenging task to rewrite the PDEs obtained this way in 
a conservative form. Eor this purpose, we employ symbolic 
integration methods that can deal with unspecified functions 
in several variables. Our approach allows us to deduce the 
mean-field equations for a general class of transportation 
processes in multiple space dimension, including the 
dynamics of multiple species that may have different size or 
shape. We illustrate our approach for a minimal model of 
pedestrian dynamics, which includes cohesion and aversion 
in bidirectional pedestrian flows. 


We have produced a prototype implementation in Math- 
ematica of the methods described in this paper, which is 
available at http://www.koutschan.de/data/meanfield/ together 
with a demo notebook. 

II. Eormal derivation of a mean-field PDE model 

FOR BIDIRECTIONAL PEDESTRIAN DYNAMICS 

We Start with a specific example to illustrate the derivation 
of a mean-field model from a discrete lattice-based approach 
in the case of bidirectional pedestrian flows. We consider two 
groups of individuals — one moving to the right, the other to 
the left. The dynamics of each individual are determined by 


cohesion and aversion — this means they try to follow and 
stay close to individuals moving in the same direction and to 
step aside when being approached by an individual moving in 
the opposite direction. We expect that these minimal dynamics 
will lead to the formation of directional lanes, a phenomenon 
that has been observed in crowded corridors, pedestrian walks 
or experiments. 


A. The microscopic model 

We start with the underlying microscopic model, i.e. a 
lattice-based approach in which we consider, for the sake of 
simplicity, a rectangle 12 C such as a corridor, partitioned 
into a square lattice of grid size h. This can be generalized 
to higher dimensions. Each lattice site {xi,yj) = {ih,jh), 
i = 0,... ,N and j = 0,... ,M can be occupied by 
an individual. We consider two groups moving in opposite 
direction — one to the right (called the reds) and one to the 
left (called the blues). The probability to find a red individual 
at time t at location {xi,yj) is given by: 

rij{t) = P(red individual is at position {xi,yj) at time t), 

where P denotes the probability. The probability for the blue 
individuals is defined analogously. We denote by 
the rate at which an individual of color c moves from (xj, yj ) 
to {xk,yi)- The transition rates for the red and blue individuals 
respectively are given by: 

7;{m}-E+L.} = (1 _ + ari+2,j), 

= (1 - Pi, j-i)ho + 71 (la) 

7-E.il-lLi+i} = (1 _ p. 
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where we write pij = ri,j + bij, and with 0 < 70 , 71,72 < 1 , 
0 < a < The prefactor {1 — p) in all terms of (1) corre¬ 
sponds to the so-called size exclusion, i.e. an individual cannot 
jump into the neighboring cell if it is occupied. We assume 
that the transition rates only depend on lattice sites in direction 
of movement, a reasonable assumption when modeling the 
movement of pedestrians. The second factors in (1) correspond 
to cohesion and aversion. Cohesion is modelled in the first line 
in each case, by introducing a factor a > 0 which increases the 
probability to move in the walking direction if the individual in 
front, i.e. in (la) at position {xi+ 2 ,yj), is moving in the same 
direction. The second and third line in each case account for 
aversion via sidestepping. If an individual, i.e. in (la) a blue 
particle located at {xi+i,yj), is approaching, the red particle 
jumps up or down (with rates 71 resp. 72 ). If 71 > 72 , the 
preference is to jump to the right with respect to the direction 
of movement, if 71 < 72 , to the left. The parameter 70 > 0 
corresponds to diffusion in the j/-direction. 


Then the evolution of the red particles is given by the so- 
called master equation 

ri,j{tk+i) = ri,j{tk) -f 

+ ( 2 a) 



Hence the probability to find a red particle at location 
{xi , Uj ) corresponds to the probability that a particle located at 
{xi-i,yj) jumps forwards (first term), particles located above 
or below, i.e. at {xi,yj±i) jump up or down (second line), 
minus the probability that a particle located at {xi,yj) moves 
forward or steps aside (third line). The evolution of the blue 
particles can be formulated analogously: 
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B. Derivation of the macroscopic model 

In the next step we formally derive the limiting mean-field 
equations as the grid size h and the time steps At tend to 
zero. Hence we consider the formal hyperbolic limit as ft, = 
At = Ax = Ay -5- 0 in Equations (2). We first substitute the 
transition rates (1) in Equation (2a) and obtain 

ri^j{tk+i)-rij{tk) = (1 - bij - rij){l + ari+ij)ri-ij 

+ (7o + 7i^i-i-i.j-i-i)(l - bi,j - rij)rij+i 
+ (7o + 72&i-i-i7-i)(l - bij - rij)rij-i 
-((1 - bi+ij - ri+ij)(l + ari+2,j) (3) 

+ (7o + 7i^i-i-i7)(l - bi,j-i - Tij-i) 

+ (7o + 72&i-ri7)(l - bi,j+i - rij+i))rij, 

and a similar equation for the evolution of the blue particles. 
Next we employ Taylor expansions up to second order of all 
the occurring probabilities. Eor example, the probability to find 
a red particle at location (xi+i,t/j) can be expanded as 

Ti+ij = Tij + hd^rij + \h‘^dlrij + 0{h^). (4) 

After expanding all probability densities we keep the terms up 
to second order and consider the formal limit as At = Ax = 
Ay 0. This leads to the following system of PDEs for the 
densities of the red and blue particles, which can be either 
obtained by tedious hand calculations, or by the computer- 
algebra methods described in Section III: 

dtr = -dx ((1 - p){l + ar)r) + (71 - 72 )^^ ((1 - p)br) 

- ^ - p){l + ar)) - 2d^{{l - p)d^r)] 

+ ^ [(71 + l 2 )dy ((1 - p)dy{rb) + brdyp) (5a) 

+ “^lody ((1 - p)dyr + rdyp) 

+ 2(71 - -f 2 )dy ((1 - p)rdxb )], 


dtb = ((1 - p)(l -f ab)b) - (71 - -f2)dy ((1 - p)br) 

- ^ - p)(l + cxb)) - 2a^((l - p)d^b)] 

+ ^ [(71 + l 2 )dy ((1 - p)dy{rb) + brdyp) (5b) 
-I- 27o(9y ((1 - p)dyb + bdyp) 

+2(71 - l2)dy ((1 - p)bd^r)]. 

The first terms on the right-hand side of (5 a) and (5b) result 
from the first-order terms in the Taylor expansion. They 
correspond to the movement of the reds and blues to the right 
and left respectively as well as the preference of either stepping 


to the right or left (depending on the difference 71 — 72). 
The second line corresponds to the second-order terms in x- 
direction, the last three lines to the second-order terms due to 
side-stepping. Note that Equations (5a) and (5b) can be written 
in a conservative form, i.e. dtr = —V • Fr and dtb = —V ■ Ft, 
for some matrices Fr and Fb. These so-called continuity 
equations are always useful as they describe the transport of a 
conserved quantity, in our case mass conservation. 

III. Algorithmic derivation of the 
MEAN-FIELD PDES 

Symbolic computation, the field of mathematics that is 
concerned with computer-implemented exact manipulation 
of mathematical expressions involving variables/symbols, is 
meanwhile a well-established area of research and has nu¬ 
merous applications. Unfortunately, it is not as widely known 
as it should be. One reason may be that some applications 
are not straightforward and require at least some insight or 
programming skills. But to those who get moderately familiar 
with symbolic computation software, it becomes an indis¬ 
pensable tool. There are plenty of general-purpose computer 
algebra systems available, the most well-known being probably 
Mathematica, Maple, and Sage. Eor our implementation we 
have chosen Mathematica. 

In this section we demonstrate how the transition from 
the discrete lattice-based model to a macroscopic PDE-based 
formulation is achieved using techniques from symbolic com¬ 
putation. 

A. Expansion 

Recall that the lattice sites are given by (xj, yj) = {ih,jh) 
for i,j € Z; in the limit ft -5- 0 one obtains the problem 
formulation for the macroscopic model. Let r = r{x,y) and 
b = b{x, y) denote the densities of red and blue particles in 
the macroscopic model. In order to perform the transition from 
partial difference equations for rtj and btj to partial differen¬ 
tial equations for r{x, y) and b{x, y), we employ formal Taylor 
expansions of the probabilities appearing in (3), as discussed 
in Section II, for example: 


Tj+ij = r + hdxT + \h'^dlr + ■ ■ 

. = V— 

ft! 

Fk 

(6a) 

ftij+i =6-1- hdyb F ^h'^dyb + ■ • 


(6b) 


Note that these calculations are done on a completely formal 
level. 

Although the expansions (6) are not available as a built- 
in command in Mathematica, it is a relatively simple task to 
implement them. We have made some effort to design our 
implementation as general as possible. This means that we do 
not fix the number of expansion variables (this corresponds 
to the dimension of the domain H). Moreover, we allow for 
discrete steps of any size, i.e, terms of the form rt+aj+b with 
a,b € Z can be handled as well. 

Eor our purposes it suffices to perform the Taylor expan¬ 
sions (6) on the master equation (3) up to second order. While 
this is a tedious calculation when done by hand, it is a trivial 



task for a computer algebra system. Still, when writing the 
result in expanded form, we obtain a huge expression for the 
right-hand side of (3): 

rdxb +ar^dxb - dxT+ bdxr+ 2rdxr+ 

-{■ (167 terms) - ^j 2 h^r{dyr){dldyb). 

Since h is considered to be very small, all terms involving 
h'^ or higher powers of h will be omitted (this corresponds 
to the polynomial reduction modulo h?). In our example 
Mathematica returns the following expression; 

rdxb -b ar'^dxb - d^r -b bd^r + 2rdxr -b 

-b (56 terms) — ^^ 2 hb'^dyr. 

Analogously, the master equation for the blue individuals 
yields a similar expression. These two PDEs in their expanded 
form cover approximately one page when printed. While this 
is still a bit unhandy for a human being, it is not at all a 
challenge for a computer. However, when we turn to more 
involved examples, it is worthwhile to spend a few thoughts 
on the implementation. As demonstrated above, expanding the 
equation after having inserted the Taylor series, results in the 
large expression (7), but most of its terms are deleted by 
the polynomial reduction, giving ( 8 ). In the present example 
this is not a big deal, but in other cases this large interme¬ 
diate expression turns out to be the bottleneck. It is then 
advantageous to systematically perform expansion-reduction 
steps on subexpressions; more precisely, to follow a bottom- 
up approach, starting at the leaves of the expression tree. 

B. Integration 

A common problem in symbolic computation is to bring the 
output of a computation into a form that is useful for a human 
being. The computing power that nowadays computers have 
allows to produce gigantic symbolic expressions without much 
effort. It can be much more difficult to extract the relevant 
information from such an output. In this spirit, we want to 
process the Taylor-expanded expressions such as ( 8 ) further, 
and rewrite them in a conservative, more compact form. 

One of the classical problems in symbolic computation 
is to determine the antiderivative of a given function. The 
first complete algorithm for the class of elementary functions 
was given by Risch [16], which was later extended to more 
general classes of functions, see for example [17]. Most of 
these algorithmic ideas found their ways into current computer 
algebra systems. 

In contrast to the classical integration problem, we shall 
consider cases that are more general, namely in the following 
three aspects. First, the given function a(a;) may not have an 
antiderivative in the prescribed class; in this case, it is desirable 
to decompose a{x) into an integrable part and remainder, i.e., 

a = dj + R, (9) 

where the remainder R is “as small as possible”. Second, 
the expressions we are dealing with involve unspecified func¬ 
tions, so that the input can be interpreted as a differential 
polynomial [18], [19]; for example, we would like to write 
the expression / • dxf as Third, our setting is 

multivariate, in the sense that we have several unspecified 


functions and several variables with respect to which we 
differentiate. 

The first algorithmic approach to the problem of inte¬ 
grating expressions with unspecified functions was proposed 
in [20], and independently for differential polynomials in [21]. 
This was generalized recently to integro-differential polynomi¬ 
als [22], [23], to differential fields [24], [25], and to fractions 
of differential polynomials [26], [27]. 

While current computer algebra systems are very good in 
computing the antiderivative of an expression involving un¬ 
specified functions (provided that it exists), the decomposition 
into an integrable part and remainder is a more delicate task. 
For example, both Mathematica and Maple correctly compute 

J{f{dx9) - ‘^idxf)^g - 2f{dlf)g) dx = 

f{dx9)-2f{dxf)g. 

In contrast, if a given expression cannot be written as the 
derivative of some other expression, then it is not at all straight¬ 
forward to obtain a decomposition of the form (9), using 
the standard integration commands provided by the computer 
algebra system. As an example, consider the decomposition 

f-dxf + f = dx{lf)+f. 

We are now going to recall the main algorithmic ideas 
how to compute a decomposition of the form (9). Fet us first 
consider a polynomial expression E in a single unspecified 
function / and its derivatives dxf,dxf ,■ ■ ■', let n denote the 
order of the highest derivative of / that appears in E. If E 
has an antiderivative, i.e., E = dxl for some polynomial 
expression I, then it is easy to see that 0”/ occurs linearly 
in E, i.e., E is quasi-linear. Hence, if 0"/ does not occur 
linearly, then the corresponding monomials are put into the 
remainder, as they cannot be integrated. Now assume that E 
is linear in 9"/. Fet m be the highest power of 9"“^/ and 
denote by u the coefficient of in E, which is 

itself a polynomial in /, dxf, • ■ •, dx~‘^f. Then integration by 
parts yields 

u-(9rV)'”(9:/) = 

Vm-fl / m d-1 

Hence the first term on the right-hand side of (10) goes into 
the integrable part, while the second term is used to replace 
u-{dx~^ f)"^{dx f) in E. After performing this step repeatedly 
(at most m times), E involves only derivatives of / up to 
order n — 1. This shows that the algorithm terminates. 

We have seen that in the case of a single unspecified 
function, there is a canonical choice which term to integrate in 
each step of the algorithm. In contrast when several unspecified 
functions are involved, the situation is less clear, as the 
following example shows: 

{dxf){dxg) = dx{f{dxg)) - f{dlg) 

= dx{{dxf)g) - {dlf)g. 

Hence one has to specify an order in which the terms are 
processed, and which at the same time doesn’t lead to infinite 



loops. The same kinds of problems are faced when the un¬ 
specified functions depend on several variables. The following 
example demonstrates the ambiguity of the decomposition in 
the case of a single unspecified function f{x,y): 

{dxf){dyf)+d^f + dyf 

= dx{f ■ dyf + f) + dy{f) - f ■ dxdyf 
= dx{f) +dy{f -dxf + f) - f ■ dxdyf. 

In our application we have to deal with several unspecified 
functions fi, ..., fk in several variables, say x,y,... ,z. So 
the question is in which order we should treat the terms of 
the input expression to obtain the desired result. One natural 
choice is to consider the variables in a fixed order as the main 
loop of the algorithm. This means that we first decompose the 
input with respect to the first variable, say dxl + R', then the 
remainder R is decomposed with respect to the next variable, 
and so on, yielding a result of the form 

dxlx + dyly + • • • + dziz + R- 

Additionally, one can also decompose I further, yielding a 
nested decomposition of the following form (we show only 
the case of a single variable): 

dx{dxi' • • {dx{I) + Rd) + ■ • • + R2) + Ri) + Ro- 

In our description of the algorithm we use the parameter d 
to specify the desired maximal integration depth of the output 
expression. 

For each integration variable, say x, we proceed as follows: 
we determine the highest derivative with respect to x that 
occurs in the input, no matter which function is involved. We 
say that the highest x-derivative is of order n if 9"/i occurs 
for some 1 < z < A:, but there is no index z such that 
for some m > 1 occurs. Then for each fi, 1 < i < k, (in 
the order as specified by the user) the terms involving d^fi 
are treated. Note that in this step derivatives with order n + 1 
can be produced, as can be seen in (11). In order to avoid that 
the algorithm runs into an infinite loop, we keep these terms, 
and continue by considering derivatives of order rz — 1. This 
algorithm is described in detail in the following pseudo-code: 


Algorithm Partialintegrate 
Input: E: differential polynomial expression 
fi, ■ ■ ■, fk- unspecified functions 
x,y,... ,z: integration variables 
d: depth 

1 : if = 0 or {x,..., 0 } = 0 or d = 0 then 

2 : return E 

3: end if 

4: i? ^ 0 

5: 1^0 

6: rz 3— highest a;-derivative that appears in E for some fi 

1 : if rz = 0 then 

8: return Partiallntegrate(ii^, {fi,---,fk),{y,---,z),d) 

9: end if 


10: for z = rz, rz — 1,..., 1 do 

11 : for j = 1,..., A: do 

12 : TO ^ HighestExponent(ii^, 9^/j) 

13: while TO > 2 do 

14: g Coefficient(E, (9 */j)’”) 

15: R^ R + g-(dlfj)"^ 

16: E ^ E - g ■ (d^f,)^ 

17: TO ^ HighestExponent(E, 

18 : end while 

19: p- 5 —Coefficient(E, 9^/j) 

20 : while p ^ 0 do 

21: TO 4— HighestExponent ( g, ^ fj ) 

22: 

23: E^E-^^ {{dlfy)g + idl-^fj)idxg)) 

24: g 4— Coefficient(E, 

25: end while 

26 : end for 

27: end for 

28: i? 4- i? -f E 

29: / 4— Partiallntegrate(/, (/i,..., fk), {x,... ,z),d — l) 
30: R ^ Partiallntegrate(i?, (/i,..., fk), (y, ■ ■ ■, z), d) 

31: return dx{I) + R 

When we apply our Mathematica implementation of algo¬ 
rithm Partialintegrate to the large expression (8) we obtain 

dt(r) = dx{r{b + r - l){ar + 1)) 

- (71 - l2)dy{br{b + r - 1 )) 

+ h(^^dx{dx{r{abr — b + ar^ — ar + 1)) -f 2rdxb) 
- (71 - l 2 )dy {r{b + r - l)dxb) 

+ li)dy{2rdyb - dy{{b - l)r)) 

+ \{li+l 2 )dy{r{ 2 b-r)dyb-dy{{b-l)br))^ 

which basically agrees with the manually derived Equa¬ 
tion (5a); recall that p = b-\-r. Comparing the two expressions 
reveals that in (5a) some remainders are not minimal according 
to algorithm Partialintegrate, but the overall expression is 
a bit more compact as it involves only factored polynomials 
inside the derivatives. 

IV. Numerical illustration for the mean-field 

MODEL 

Finally we would like to illustrate the behavior of so¬ 
lutions to (5) with numerical simulations. We consider the 
system (5) on 11 x (0, T), where in our computational examples 
n = [—Lx, Lx] X \—Ly, Ly] C ^ 2 ^ 3 , corresponds 

to a corridor. As individuals cannot penetrate the walls, we 
set no flux boundary conditions on the top and bottom. At the 
entrance and exit of the corridor, i.e at x = ELx, we assume 
periodic boundary conditions. For all numerical simulations 
we used the COMSOF Multiphysics Package with quadratic 
finite elements. We set H = [0,1] x [0,0.1], choose a mesh 
of 608 triangular elements and a BDF method with maximum 
time step 0.1 to solve the discretized system. The first example 
models system (5) in the case where we have no cohesion 
and no preference for stepping to one side, i.e. a = 0 and 





7 ■= 7 i = 72 - In the second example we consider system (5) 
with the special scaling 71—72 = 0{h) including cohesion and 
aversion. For this particular scaling the first-order hyperbolic 
terms in y-direction are of the same order as the diffusion in 
this direction while the mixed derivative terms, i.e. the terms 
which involve derivatives with respect to x and y, are of order 
0{h?) and can be neglected. In both simulations we start 
with a perturbed configuration of the steady state, i.e. constant 
densities for r and b, and study if the densities return to the 
constant steady state or to another more complex stationary 
configuration. 

1) Example 1: Let 70 = 0.1, 7 = 0.2 and h = 0.3. As 
initial values we choose small perturbations of an equilibrium 
state, i.e. 

ro(a;, y) = 0.4 + 0 . 02 sin( 7 rx) cos 

'" 0 - 1 ^ ( 12 ) 

ba{x,y) = 0.4 — 0.02sin( 7 rai) cos ■ 

The initial value tq and the solution tt to the system (5) 
at time T = 5 is visualized in Figure 1. The corresponding 
density of blue individuals show the same behavior, i.e. the 
densities return to the constant equilibrium solution. In this 


Time=0 s Surface: Dependent variable r (m) 



(a) Initial distribution of reds 

Time=5 s Surface: Dependent variable r (m) 



(b) Distribution of reds at time T = 5 


Fig. 1. Solution to system (5) with no cohesion and aversion. 

example we do not observe the formation of directional lanes 
as the solutions return to the equilibrium state quickly. 

2) Example 2: Let 70 = 0.001, 71 = 0.5, 72 = 0.4, 
a = 0.2 and h = 0.1, i.e. 71 —72 = 0{h). As initial values we 
choose again ( 12 ), i.e. small perturbations of an equilibrium 
state. Figure 2 shows the solution and bx to system (5) 
at time T = 5. In this example we observe lane formation. 



(a) Density of reds at time T = 5 



(b) Density of blues at time T = 5 
Fig. 2. Solution to system (5) with cohesion and aversion. 

Since 71 > 72 individuals have a tendency to step to the 
right. This tendency can also be observed in the formation 
of the directional lanes. The red individuals concentrate on 
the bottom of the corridor, whereas the blue individuals move 
to the top. 

V. Conclusion and outlook to further 

APPLICATIONS 

In this paper, we have presented a symbolic approach 
to derive the mean-field PDFs from lattice-based models. 
We demonstrated the methods in terms of one example in 
pedestrian dynamics, but the algorithm may also be applied to 
other examples. In [1] different motility mechanisms on regular 
lattices are introduced, which result in nonlinear diffusion 
equations with different diffusivities. The authors considered 
various motilities based on attraction or repulsion, i.e. the 
transition rate to move away from a neighbouring individual 
increases or decreases respectively. For example, in a minimal 
model the transition rate is given by 

= (l-c,+i)(l-ac,_i), 

where Ci again denotes the probability that the lattice site Xi 
is occupied. Hence the transition to move from Xi to Xi+i 
is reduced if the neighbouring site Xi-i is occupied. This 
phenomenon is known as adhesion and results in a nonlinear 
diffusion model for the cell density c = c{x, t) of the form 

dtc = dx{D{c)dxc), 

with a diffusivity of the form D(c) = 3a(c — |)^ -I- 1 — |a, 
see [28]. Again in [1] several other transition rates were 





























proposed, which lead to different nonlinear diffusitivies, see 
Table 1 there. Using our implementation, the entries in that 
table can be generated automatically. For example, we have 
tried one of their most complicated models [1, Equation (13)], 
which combines contact-forming or contact-breaking interac¬ 
tions with contact-maintaining interactions. Our implementa¬ 
tion correctly derives the diffusivity given in [1, Equation (14)], 
where we have chosen the two-dimensional square lattice with 
Moore interacting neighborhoods. 
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